%% First load data arrays for analysis
clear;
fnames = dir("monod*.mat");
for k=1:length(fnames)
    load(fnames(k).name);
end

%% Define molar mass
molmass_fructose = 180; % Molar mass of fructose
clc;
close all;
%% Do 37C
tempmat = monod_mat_fructose_25C;
smat = size(tempmat);

conc1 = tempmat(:,1); % Convert to mM
grmax1 = nanmean(tempmat(:,2:3),2);
grmax_err1 = nanstd(tempmat(:,2:3),0,2);

% Do fit on corrected concentration (in mM)
figure;
errorbar(conc1,grmax1,grmax_err1, 'bo', 'linewidth', 1);
    
hold on;

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')
box off;

% Monod fit
[f1,f2] = monod_fit_weighted(conc1,grmax1,grmax_err1);
ci = confint(f1);

KM = f1.b;
KMerr = abs(ci(2,2) - KM);

disp('Fit and error of KM (mM)');
disp(KM);
disp(KMerr);


%% Do 37C
tempmat = monod_mat_fructose_30C_v1;
smat = size(tempmat);

conc1 = tempmat(:,1); % Convert to mM
grmax1 = nanmean(tempmat(:,2:3),2);
grmax_err1 = nanstd(tempmat(:,2:3),0,2);

% Do fit on corrected concentration (in mM)
hold on;
errorbar(conc1,grmax1,grmax_err1, 'ko', 'linewidth', 1);
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')
box off;

% Monod fit
% Monod fit
[f1,f2] = monod_fit_weighted(conc1,grmax1,grmax_err1);
ci = confint(f1);

KM = f1.b;
KMerr = abs(ci(2,2) - KM);

disp('Fit and error of KM (mM)');
disp(KM);
disp(KMerr);

%% Do 37C
tempmat = monod_mat_fructose_37C;
smat = size(tempmat);

ind = 1:4;
conc1 = tempmat(ind,1); % Convert to mM
grmax1 = nanmean(tempmat(ind,2:3),2);
grmax_err1 = nanstd(tempmat(ind,2:3),0,2);

% Do fit on corrected concentration (in mM)
hold on;
errorbar(conc1,grmax1,grmax_err1, 'ro', 'linewidth', 1);
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
[f1,f2] = monod_fit_weighted(conc1,grmax1,grmax_err1);
ci = confint(f1);

KM = f1.b;
KMerr = abs(ci(2,2) - KM);

disp('Fit and error of KM (mM)');
disp(KM);
disp(KMerr);